February 1, 2008 4:8 WSPC/Trim Size: 9in x 6in for Proceedings nishimura 



m 
o 
o 

(N| ; FACTORIZATION METHOD FOR SIMULATING 

QCD AT FINITE DENSITY 

a ■ 



> 

in 



JUN NISHIMURA* 

Department of Physics, Nagoya University, 

Nagoya 464-8602, Japan 
E-mail: nisimura@eken.phys.nagoya-u.ac.jp 



f— *^ ■ We propose a new method for simulating QCD at finite density. The method is 

\f} ' based on a general factorization property of distribution functions of observables, 

f— *^ _ and it is therefore applicable to any system with a complex action. The so-called 

, overlap problem is completely eliminated by the use of constrained simulations, 

i We test this method in a Random Matrix Theory for finite density QCD, where 

we are able to reproduce the exact results for the quark number density. 

Qh, 1. Introduction 

■ Recently there are a lot of activities in QCD at finite density, where in- 

teresting phases such as a superconducting phase have been conjectured to 
appear . At zero chemical potential Monte Carlo simulations of lattice 
QCD enables nonperturbative studies from first principles. It is clearly de- 
sirable to extend such an approach to finite density and explore the phase 
diagram of QCD in the T(temperature)-/i(chemical potential) plane. The 
main obstacle is here that the Euclidean action becomes complex once the 
chemical potential is switched on. 

Nevertheless QCD at finite density has been studied by various ap- 
proaches with exciting conjectures. First there are perturbative studies 
which are valid in the fi — > oo limit 2,3 . Refs. [4] and [5] uses effective the- 
ories with instanton-induced four-fermi interactions. As for Monte Carlo 
studies two directions have been pursued so far. One is to modify the 
model so that the action becomes real. This includes changing the gauge 
group 6 from SU(3) to SU(2), and introducing a chemical potential with 
opposite signs for up and down quarks 7 . The other direction is to explore 



*Work partially supported by Grant-in-Aid for Scientific Research (No. 14740163) from 
the Ministry of Education, Culture, Sports, Science and Technology. 



1 



February 1, 2008 4:8 WSPC/Trim Size: 9in x 6in for Proceedings 



nishimura 



2 

the large T and small /x regime of lattice QCD, where the imaginary part 
of the action is not very large 8,9 . These studies already produced results 
relevant to heavy ion collision experiments, but more interesting physics 
will be uncovered if larger /i regime becomes accessible by simulations. 

In Ref. [10] we have proposed a new method to simulate systems with a 
complex action, which utilizes a simple factorization property of distribu- 
tion functions of observablcs. Since the property holds quite generally, the 
approach can be applied to any system with a complex action. The most 
important virtue of the method is that it eliminates the so-called overlap 
problem, which occurs in the standard re-weighting method. Ultimately we 
hope that this method will enable us, among other things, to explore the 
phase diagram of QCD at finite baryon density. As a first step we have 
tested 11 the new approach in a Random Matrix Theory for finite density 
QCD 12 , which can be regarded as a schematic model for QCD at finite 
baryon density. 

2. Random Matrix Theory for finite density QCD 

The Random Matrix Model we study is defined by the partition function 

Z = JdWe- Nt < wiw UetD , (2.1) 
where W is a N x N complex matrix, and D is a 2N x 2N matrix given by 

D=(.™ iW + fl ) . (2.2) 

The parameters m and \x correspond to the 'quark mass' and the 'chemical 
potential', respectively. In what follows we consider the massless case (m = 
0) for simplicity and we focus on the 'quark number density' defined by 

»=±*frD-*), 7.= (;j). (2-3) 

The vacuum expectation value (VEV) of the quark number density is ob- 
tained exactly by [13] and in particular in the large- TV limit 12 

lim <„) = ( |° r " < * (2.4) 
JV^oo [ 1/fj, for /x > n c , 

where /x c is the solution to the equation l+ii 2 +ln(it 2 ) = 0, and its numerical 
value is given by /i c = 0.527 • • •. We find that the quark number density 
(u) has a discontinuity at \x — \x c . Thus the schematic model reproduces 
qualitatively the first order phase transition expected to occur in 'real' QCD 
at nonzero baryon density. 
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3. The complex action problem 

Let us first rewrite (2.1) as 



dWe~ So+ir , (3.1) 

where we have introduced So and T by 

So = NtT(W*W)-ln\detD\ (3.2) 
detD = e lT \detD\ . (3.3) 

In this form it becomes manifest that the system has a complex action, 
where the problematic imaginary part T is given by the phase of the fermion 
determinant. Since the weight e~ Sa+ir in (3.1) is not positive definite, we 
cannot regard it as a probability density. Hence it seems difficult to apply 
the idea of standard Monte Carlo simulations, which reduces the problem 
of obtaining VEVs to that of taking an average over an ensemble generated 
by the probability density. 

Let us define the so-called phase quenched partition function 

Z Q = J dWe~ NtI ^ w ^ \detD\ = J dWe~ s ° . (3.4) 

Since the system (3.4) has a positive definite weight, the VEV ( • )o as- 
sociated with this partition function can be evaluated by standard Monte 
Carlo simulations. Then one can use the standard re- weighting formula 

<"> = W <3 ' 5) 

to obtain the VEV (v) in the full model (3.1). The problem with this 
method is that the fluctuations of the phase V in (3.5) grows linearly with 
the size of the matrix D, which is of O(N). Due to huge cancellations, both 
the denominator and the numerator of the r.h.s. of (3.5) vanish as e _const - N 
as N increases, while the 'observables' e* r and ue lT are of 0(1) for each 
configuration. As a result, the number of configurations required to obtain 
the VEVs with some fixed accuracy grows as e constJv . 

In fact we may simplify the expression (3.5) slightly by using a symme- 
try. We note that the fermion determinant det D, as well as the observable 
v, becomes complex conjugate under the transformation 



W^-W , 



(3.6) 
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while the Gaussian action remains invariant. From this we find that 

M = K) + i (ia) (3.7) 
_ (^R.cosr) fosinr) r „ ^ 

(cosr) (cosr) 

where i/r and z^i denote the real part and the imaginary part of v, respec- 
tively. This simplification, however, does not solve the problem at all, since 
cos T and sin T flip their sign violently as a function of the configuration W . 
Note that both terms in the r.h.s. of (3.7) are real, meaning in particular 
that their sum (y) is also real. 

The model (3.4) is solvable in the large- N limit 12 and one obtains 

lnnMo = ((\ 0.9) 
a^oo' ( 1/fj, for n > 1 . 

In this case the VEV of the quark number density is a continuous function 
of the chemical potential p, unlike in (2.4). Thus the first order phase 
transition in the full model (3.1) occurs precisely due to the imaginary part 
r of the action. Note also that the symmetry under (3.6) implies 

(^i)o = ; <i/r>o = Mo ■ (3.10) 
4. The factorization method 

In this section, we explain how the factorization method 10 can be used to 
obtain the VEVs (z/r) and (i^). The fundamental objects of the method 
are the distribution functions 

Pi (x) d ^ (S(x - i/O) (4.1) 

p\ 0) (x)^{5{x- Vi )) i = R,I (4.2) 

defined for the full model and for the phase quenched model respectively. 
The important property of these functions is that they factorize as 

Pi(x) = ^pf\x)<pi{x) i = R,I, (4.3) 

where the constant C is given by C = (e iT ) . The 'weight factor' ipi(x) 
represents the effect of T, and it can be written as a VEV 

Viix)^ 1 (e ir ) itX (4.4) 
with respect to a yet another partition function 

Zi{x) = J dWe~ s « S(x - Vi) . (4.5) 
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The <5-function represents a constraint on the system. In actual simulation 
we replace the (5-function by a sharply peaked potential. We refer the reader 
to Ref. [11] for the details. 

Using the symmetry under (3.6), the formulae for (v) is nothing but 
(3.8), where (^R,cosr) , (^isinr} and (cosr) are replaced by 

/oo 
dxxp^\x) w R (x) , (4.6) 
-OO 

(i/isinr) = 2 / dx x p[ 0) {x) wi(x) , (4.7) 
Jo 

/oo 
dxp { °\x)w R (x) . (4.8) 
-oo 

The weight factors Wi(x) are defined by 

w R (x) = (cosT) R , x ; wi{x) d = (smT)i. x . (4.9) 

One of the virtues of the method can be seen from (4.6)~(4.8). If we are 
to obtain the VEVs on the l.h.s. by directly simulating the system (3.4), for 
most of the time we sample configurations whose Vi takes a value close to the 
peak of p\°\x). However, from the r.h.s. of the formulae, it is clear that we 
have to sample configurations whose vi takes a value where p] (x)\wi(x)\ 
becomes large, in order to obtain the VEVs accurately. In general these 
two regions of configuration space have little overlap, which becomes expo- 
nentially small as the system size increases. The present method resolves 
this 'overlap problem' completely by 'forcing' the simulation to sample the 
important region. 



Table 1. Results of the analysis of (u) described in the text. Sta- 
tistical errors computed by the jackknifc method arc shown. The last 
column represents the exact result for (u) at each /i and N . For /i = 0.2 
the exact result is (u) = —0.2 with an accuracy better than 1 part in 
10" 9 . 





N 


<^R> 


i(vi) 


<") 


(u) (exact) 


0.2 


8 


0.0056(6) 


-0.1970(5) 


-0.1915(7) 


-0.20000. . . 


0.2 


16 


0.0060(4) 


-0.1905(13) 


-0.1845(13) 


-0.20000. . . 


0.2 


24 


0.0076(9) 


-0.1972(14) 


-0.1896(17) 


-0.20000. . . 


0.2 


32 


0.0021(8) 


-0.1947(19) 


-0.1927(25) 


-0.20000. . . 


0.2 


48 


0.0086(37) 


-0.2086(54) 


-0.2000(88) 


-0.20000. . . 


1.0 


8 


0.8617(10) 


0.1981(13) 


1.0598(12) 


1.066501. . . 


1.0 


16 


0.8936(2) 


0.1353(6) 


1.0289(5) 


1.032240. . . 


1.0 


32 


0.9207(1) 


0.0945(2) 


1.0152(3) 


1.015871. . . 
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5. Reproducing exact results by the new method 

In Table 1 we show our results for two values of p, p = 0.2 and p = 1.0, 
which are on opposite sides of the first order phase transition point p = 
p c = 0.527 • • •. They are in good agreement with the exact results, and the 
achieved values of N are large enough to extract the large ./V limit. 

Note that (i/r) <~ at p = 0.2. Thus the main contribution to (y) 
comes from the imaginary part (vi), which is in sharp contrast to the results 
(3.10) for the phase quenched system. This result comes about because the 
sign change of u>r(x) occurs near the peak of p^\x), so that the product 
\x)wr(x) has a positive regime and a negative regime, which cancel 
each other in (4.6). For p = 1.0, on the other hand, u>r(x) is approximately 
constant in the region where p^ (x) is peaked, so the shape of the product 
P^\x)wr(x) is similar to p^(x). The main contribution to (y) comes from 
the real part (vr), and moreover, it is close to (^r)o- 
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Figure 1. The weight factor wr (x) is plotted against x for TV 
behavior changes drastically as /i crosses the critical point. 



at various [i. The 



In Fig. 1 we plot wr(x) for N — 8 at various p. It is interesting that the 
wr(x) changes from positive to negative for /i < p c , but it changes from 
negative to positive for p > p c . (Similarly w\(x) is positive at x > for 
p < p c , but it is negative at x > for p > p c .) Thus the behavior of Wi(x) 
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changes drastically as the chemical potential fi crosses its critical value /i c . 
These results provide a clear understanding of how the first order phase 
transition occurs due to the effects of T. 

6. Applications to other systems with complex actions 

The method 14 proposed for simulating ^-vacuum like systems can be re- 
garded as a special case of the factorization method. A simplified version 
of the method was sufficient because the observable was identical to the 
imaginary part of the action. The essence of the factorization method is 
that it avoids the overlap problem by the use of constrained simulations. 
In Ref. [14] promising results for 2d CP 3 are also reported. 

In Ref. [10] the method has been used to study the dynamical generation 
of space time in superstring theory based on its matrix model formulation 
15 . There the method becomes even more powerful since the distribution 
functions turn out to be positive definite. In this case the scaling property 
of the weight factor enables extrapolations to larger system size. 

We hope that the factorization method is useful also for studying other 
interesting systems with complex actions such as Chern-Simons theories, 
chiral gauge theories, strongly coupled electron systems etc. 
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